Draft version March 5, 2013 

Preprint typeset using KTpX style emulateapj v. 08/13/06 



FORMATION OF TURBULENT AND MAGNETIZED MOLECULAR CLOUDS VIA ACCRETION FLOWS OF HI 

CLOUDS 

TSUYOSHI INOUE 1 AND SHU-ICHIRO INUTSUKA 2 
Draft version March 5, 2013 

ABSTRACT 

Using three-dimensional magnetohydrodynamic simulations including the effects of radiative cool- 
ing/heating, chemical reactions, and thermal conduction, we investigate the formation of molecular clouds 
in the multi-phase interstellar medium. We consider the formation of molecular clouds due to accretion of Hi 
clouds as suggested by recent observations. Our simulations show that the initial Hi medium is compressed and 
piled up behind the shock waves induced by accretion flows. Since the initial medium is highly inhomogeneous 
as a consequence of the thermal instability, a newly-formed formed molecular cloud becomes very turbulent 
owing to the development of the Richtmyer-Meshkov instability. The structure of the postshock region is com- 
posed of dense cold gas (T < 100 K) and diffuse warm gas (T > 1, 000 K), which are spatially well mixed 
due to turbulence. Because the energy source of the turbulence is the accretion flows, the turbulence becomes 
highly anisotropic biased towards the direction of the accretion flows. The kinetic energy of the turbulence 
dominates the thermal, magnetic, and gravitational energies throughout the entire 10 Myr evolution. However, 
the kinetic energy measured using CO-fraction-weighted densities is comparable to the other energies at t ~ 5 
Myr, once the CO molecules are sufficiently formed owing to UV shielding. This suggests that the true kinetic 
energy of turbulence in molecular clouds as a whole can be much larger than the kinetic energy of turbulence 
estimated using line-widths of molecular emission. We find that clumps in a molecular cloud show a statisti- 
cally homogeneous evolution as follows: the typical plasma (i of the clumps is roughly constant ((3) ~ 0.4; 
the size-velocity dispersion relation is At> ~ 1.5 km s _1 (7/1 pc) 5 , irrespective of the density; the clumps 
evolve toward magnetically supercritical, gravitationally unstable cores; the clumps seem to evolve into cores 
that satisfy the condition for fragmentation into binaries. These statistical properties may represent the initial 
conditions of star formation. 

Subject headings: galaxies: ISM — ISM: clouds — ISM: magnetic fields — ISM: molecules — stars: forma- 
tion 



1. INTRODUCTION 

It is well known that molecular clouds are the sites of 
present-day star formation. However, our understanding of 
the physical conditions of molecular clouds is very limited. 
Investigations of the formation of molecular clouds should be 
one of the most promising way to reveal the physical condi- 
tions of molecular clouds. Recent progress in numerical sim- 
ulations of the interstellar medium has revealed many vital as- 
pects of interstellar cloud formation. Using one-dimensional 
hydrodynamics simulations, Hennebelle & Perault (1999) and 
Koyama & Inutsuka (2000) showed that shocked warm neu- 
tral medium (WNM or diffuse intercloud medium; T <~ 10 4 
K) evolves into cold neutral medium (CNM or interstellar 
cloud; T < 10 2 K) owing to radiative cooling. In the follow- 
up two-dimensional simulation, Koyama & Inutsuka (2002) 
found that the thermal instability (Field 1965, Field et al. 
1969) that is caused by runaway cooling in shocked WNM 
generates highly clumpy small-scale clouds with velocity dis- 
persion that are supersonic with respect to the CNM. 

Since supersonic turbulence is a universal feature of in- 
terstellar clouds (Larson 1981, Heiles & Troland 2005), the 
formation of turbulent clouds without ad-hoc external forc- 
ing is very advantageous to the study of the physical nature 
of clouds. Hennebelle & Audit (2007) and Hennebelle et 
al. (2007) showed that Hi clouds formed by thermal insta- 
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bility can explain various observational characteristics. In ad- 
dition to the turbulence, magnetic fields are also an impor- 
tant ingredient of the interstellar clouds (Heiles & Crutcher 
2005). Using magnetohydrodynamic (MHD) simulations, In- 
oue & Inutsuka (2008, 2009) showed that molecular clouds 
can be formed behind a shock wave if the shock propagates 
almost parallel to the mean magnetic field, otherwise only Hi 
clouds are formed since magnetic pressure prevents the ac- 
cumulation of the gas (see also Hennebelle & Perault 2000, 
Hartmann et al. 2001, Inoue et al. 2007). Recently, em- 
ploying three-dimensional MHD simulations, Hennebelle et 
al. (2008), Heitsch et al. (2009), Banerjee et al. (2009), 
and Vazquez-Semadeni et al.(2011) studied the formation of 
molecular clouds due to converging flows of WNM along 
magnetic field lines. Note that in the above-mentioned studies 
the gas is assumed to be atomic, and optically thin radiative 
cooling/heating are employed in the simulations. Thus, the 
effects of chemical evolution including the effects of ultravi- 
olet (UV) shielding and the resulting modifications to radia- 
tive cooling/heating mechanisms should be taken into account 
(e.g. Glover et al. 2007, 2010). 

Recent observations of nearby galaxies, molecular clouds 
are suggested to be formed from Hi clouds with n <~ 10 cm~ 3 
(Blitz et al. 2006, Fukui et al. 2009). A detailed statisti- 
cal study of molecular clouds in the Large Magellanic Cloud 
(LMC) showed that the timescale of molecular cloud evolu- 
tion due to Hi cloud accretion is approximately <~ 10 Myr 
(Fukui et al. 2008, Kawamura et al. 2009). As discussed 
by Hartmann et al. (2001) and Inoue & Inutsuka (2009), 
since the ISM can only be accumulated efficiently parallel 
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to the mean magnetic field, forming a cloud with visual ex- 
tinction (A v ) > 1 (N H > 2 x 10 21 cirT 3 ) by piling up 
the WNM (n ~ 1 cm~ 3 ) takes t {oIm > 80 Myr (A v /1) 
(wflow/8kms _1 ) _1 (n/lcm -3 ) -1 , even if the WNM is ac- 
cumulated at the sound speed Vfl ow = c s,wnm ~ 8 km s _1 . 
This is an order of magnitude larger than the observational 
estimation of ~ 10 Myr. Thus, in this paper, we examine the 
formation of molecular clouds by the accretion of Hi clouds 
denser than the WNM, as suggested by observations. 

The plan of this paper is as follows: in §2, we provide a 
possible scenario for molecular cloud formation and the nu- 
merical setup of our simulations; the results of the simulations 
are shown in §3; finally, in §4, we summarize our findings and 
discuss their implications. 

2. FORMATION SCENARIO AND NUMERICAL SETUP 
2.1. Scenario 

In Inoue & Inutsuka (2008, 2009), we showed for the case 
in which the angle between the shock normal and the mean 
magnetic field is not parallel that the WNM compressed by 
an interstellar shock wave (such as a supernova blast wave, 
super-bubble, or galactic spiral shock) is cooled roughly iso- 
chorically and generates thermally unstable gas. The resulting 
thermally unstable medium evolves into a complex of sheet- 
like Hi clouds with n ~ 30 cm -3 embedded in the diffuse 
WNM. Since the total mass of Hi clouds formed behind a 
super-bubble can be huge, they could be a mass reservoir of 
the building blocks of molecular clouds. According to In- 
oue & Inutsuka (2009), the average number density of the Hi 
medium composed of Hi clouds and WNM generated behind 
a shock is ~ 10 cm~ 3 , which agrees well with the observed 
Hi mass reservoir of molecular clouds (Fukui et al. 2009, 
Dawson et al. 201 1). 

If the Hi reservoir is swept up along the mean magnetic field 
by shock waves, e.g. due to a galactic spiral shock or super- 
bubble, molecular clouds can be formed. In this case, because 
the mean density of the medium is an order of magnitude 
larger than the WNM, the timescale of mass accumulation to 
form a cloud of (A v ) > 1 becomes tf oxm > 8 Myr (A v /1) 
(fflow/8kms _1 ) _1 (rt/10cm -3 ) -1 , consistent with the evo- 
lutionary timescale of observed molecular clouds (Kawamura 
et al. 2009). In addition, if the Hi reservoir is massive 
enough to be gravitationally unstable, the formation of molec- 
ular clouds due to gravitational contraction along mean mag- 
netic field lines would also be possible. Thus, in this paper, 
we examine the formation of molecular cloud due to the tran- 
sonic accumulation of Hi clouds along magnetic field lines. 

2.2. Formation of Hi Reservoir 

In this section, we provide the computational setup for 
forming the Hi reservoir for molecular cloud formation. The 
detailed modeling of molecular cloud formation from this Hi 
reservoir is presented in the next section. The basic equations 
to be solved are the ideal MHD equations including the effects 
of chemical reactions, radiative cooling/heating, and thermal 
conduction. 
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where the subscript i denotes the chemical species p, H, H 2 , 
He, He + , C, C + , and CO, /; is the chemical reaction term 
for species i, and L is the net cooling rate per unit volume. 
The other symbols have their usual meaning. The chemical 
reaction coefficients and the net cooling rate depend on tem- 
perature (T), the strength of the external radiation field (Go), 
and the column densities of chemical species and dust (iVj)- 
In this paper, we consider a requisite minimum set of chemi- 
cal reactions and cooling/heating processes in order to follow 
the formation of an opaque molecular cloud from an optically 
thin Hi medium. We summarize the chemical reactions and 
cooling/heating processes in Tables 1 and 2, respectively. 

The basic MHD equations are solved using a second-order 
Godunov-type finite volume scheme (van Leer 1979) that em- 
ploys an approximate Riemann solver developed by Sano et 
al. (1999). The consistent method of characteristics with con- 
strained transport is used for solving the induction equation 
(Clarke 1996, see also Evans & Hawley 1988, Stone & Nor- 
man 1992). The integration of the ideal MHD part of the ba- 
sic equations is done in a conservative fashion. The cooling, 
heating, and thermal conduction terms are solved using the 
second-order explicit method. We use the piecewise exact so- 
lution method developed by Inoue & Inutsuka (2008), an un- 
conditionally stable second-order method originally adopted 
for solving friction terms in two-fluid MHD equations, to cal- 
culate the chemical reaction terms. We confirmed the accu- 
racy of this method by solving the chemical reaction equa- 
tions in one-zone test problems. We find that the piecewise 
exact solution method is as accurate as the Crank-Nicolson 
scheme, even though it does not require iterative procedures 
or matrix inversion (Kamata 2012). The timestep for the nu- 
merical integration is determined by the minimum timestep 
due to the CFL condition for the MHD fast mode, one-fifth 
of the local cooling time, and the stability condition for the 
parabolic conduction equation. 

We prepare an Hi medium composed of Hi clouds and 
WNM as building blocks of a molecular cloud. In this stage, 
we assume optically thin cooling, heating and chemical reac- 
tions, i.e., the effects of UV shielding in the chemical reac- 
tion and heating/cooling terms are switched off. The back- 
ground UV field strength is set to the Habing flux Go = 1 
(= 1.6x 10~ 3 erg cm~ 2 s _1 ). The thermal and chemical equi- 
librium state of this optically thin medium is shown as a solid 
line in the number density-pressure plane in panel (a) of Fig- 
ure [T] We setup an initial state of thermally unstable gas with 
density fluctuations, the spectrum of which is a power law 
with the Kolmogorov spectral index (P3D oc fc~ n / 3 ). Ther- 
mal pressure and average number density are set to p/k-Q = 

5.2 x 10 3 K cm~ 3 and (n) = J2( n i) = 5-2 cm~ 3 , where 
the initial abundances are ih = n n/J2 n i — 0.91, x p = 
1.7 x 10~ 3 , x H , = 1-3 x 10~ 6 , xhc = 0.090, x Hc + = 

1.3 x 10~ 4 , x c l = 1.4 x 10~ 4 , x c = 6.0 x 10~ 9 , and 
^co = 5.7 x 10~ 18 . The initial mean molar weight and the 
total mass in the numerical domain are m mcan = 1.27 m p and 
A/tot = 1305 M so \, respectively. The mean density and pres- 
sure of this initial state is shown as a cross in Figure |TJa). We 
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TABLE 1 
Chemical Reactions 





Note 


R 1Y 1 vt^ n p 

IVC1C1C1H-C 


Cosmic ray ionization 


Ch = 3.0 x 10 17 s-\CHe = 1-087 Ch.Cc = 3.84( 


Kh 1 


H2 formation on grains 


dust temperature T c i = 10 K is used 


2 


H2 photo-dissociation 


depends on UV strength Go, column density Afjj, , 






and visual extinction Av 


3 


H2 and CO dissociation 


collisions with e, p, and H 


4 


H2 and CO dissociation 


recombination of He+ 


1 


H, He, and C ionization 


collisions with e, p, H, and H2 


1.4 


H+, He" 1 ", and C + recombination 




4,5 


CO formation 


depends on Go and Av 


6 


CO photo-dissociation 


depends on Go, Nn 2 , Nco an d Av 


6,7 


C photo-ionization 


depends on Go, Ac, ^Vh 2 an d Av 


2 


References. — (1) Millar et al. 


1997; (2) Tilelens & Hollenbach 1985; (3) Draine & Bertoldi 1996; (4) 


Hollenbach & McKee 1989; (5) Shapiro & Kang 1987; (6) Nelson & Langer 1997; (7) Lee et al. 1996 




TABLE 2 




Cooling and Heating Processes 




Process 


Note 


Reference 


Photoelectric heating by PAHs 


depends on Go and Av 


1,2 


Cosmic ray heating 




3 


H2 photo-dissociative heating 


depends on Go, Ah 2 , and Av 


4 


Ly-ct cooling 




5 


C+ cooling (158/^m) 


level population and escape probability 






as a function of N c + are taken into account 


6 


O cooling (63/im) 


escape probability as a function of Nq = xq Ntot 






is taken into account (xq = 3.2 X 10 — 4 ). 


2.6 


CO ro-vibrational cooling 


escape probability as a function of A^co 






is taken into account 


7,8,9 


Cooling due to rec. of 






electrons with grains & PAHs 




1 


Thermal conduction 


collisions of H and H2 


10 


REFERENCES. — (1) Bakes & Tielens 1994; (2) Wolfire et al. 2003; (3) Goldsmith & Langer 


1978; (4) Black & Dalgano 1977; (5) Spitzer 1978; (6) de Jong et al. 1980; (7) Hollenbach & 


MacKee 1979; (8) Hosokawa & 


Inutsuka 2006; (9) Hollenbach & McKee 1989; (10) Parker 1969 



it to evolve towards the WNM phase. Linear stability analyses 
of the thermal transition layer between an Hi cloud and WNM 
have shown that the transition layer is always unstable with re- 
spect to corrugational perturbations (Inoue et al. 2006, Stone 
& Zweibel 2009). Thus, the resultant Hi clouds are highly 
corrugated. The sheet-like morphology of the Hi clouds, their 
magnetization (~ 5 /iG), typical density (n ~ 50 cm -3 ), and 
temperature (T ~ 100 K) agree well with observed Hi clouds 
in the ISM (Heiles & Troland 2003). In addition, the mean 
density of the Hi reservoir composed of Hi clouds and the 
WNM is also consistent with observed Hi reservoirs (Fukui 
et al. 2009). As we will show in §3, the mass spectrum of the 
Hi clouds is fitted by a power-law of dN{M)/dM oc M -1 ' 8 
that shows good agreement with the theoretical mass spec- 
trum of Hi clouds formed as a result of the thermal instability 
(Hennebelle & Audit 2007, Hennebelle & Chabrier 2008). 

2.3. Molecular Cloud Formation 

In this section, we provide our numerical setup for molec- 
ular cloud formation. We use the Hi medium obtained in the 
previous section as an initial condition, and consider the con- 
verging flows of the Hi reservoir along the mean magnetic 
field (±2>direction). For this purpose, we set the initial ve- 
locity to v x = u cnv tanb-(— x/1 pc), where the center of the 
numerical domain is defined as x = y = z = 0. In this 
paper, we examine the case of a converging flow velocity of 



use a cubic numerical domain of (20 pc) 3 in volume which 
we divide into 1024 3 uniform cells, such that the numerical 
resolution is ~ 0.02 pc. The initial medium is uniformly 
magnetized by a field of strength 5.0 /iG oriented in the +x- 
direction. As shown by Inoue & Inutsuka (2009), such a ther- 
mally unstable medium can be formed behind a shock wave 
that propagates in turbulent WNM. Periodic boundary condi- 
tions are used for this stage of Hi cloud formation. 

The three-dimensional structure of the density resulting 
from 8 Myr of integration (several cooling times for the ini- 
tial thermally unstable gas) is shown in the top panel of Fig- 
ure]^ Regions in blue are Hi clouds with n > 10 cm~ 3 and 
T ~ 100 K that are formed as a consequence of the thermal 
instability. Magnetic field lines are shown as black lines. A 
two-dimensional slice is also shown, in which regions in red 
represent n < 3 cm -3 , regions in yellow and green repre- 
sent 3 cm -3 < n < 10 cm— 3, and regions in blue represent 
n > 10 cm~ 3 . The probability distribution function (PDF) of 
the number density is plotted in the bottom panel of Figure [2] 
Condensation driven by the thermal instability arises along 
the magnetic field lines, since condensation perpendicular to 
the field is suppressed due to the enhancement of magnetic 
pressure (Hennebelle & Perault 2000, Inoue et al. 2007, In- 
oue & Inutsuka 2008, 2009). This results in the formation of 
sheet-like clouds. Note that the gas in the intercloud medium 
is also thermally unstable due to runaway heating that causes 
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FIG. 1 . — Panel (a): Thermal and chemical equilibrium curve for media of 
Av=0 (optically thin, solid) and Av=l (dashed). Dotted lines are isotherms 
of T = 10 1 , 10 2 , 10 3 , and 10 4 K. Panel (b): Cooling and heating rates of 
the thermal and chemical equilibrium states for Av=0 (optically thin, solid) 
and Av=l media (dashed). The symbol Ly-a denotes Ly-a cooling, CII de- 
notes Cn fine structure line (158/^m) cooling, OI denotes Oj fine structure 
line (63/xm) cooling, DREC denotes the effect of cooling due to the recom- 
bination of electrons with grains and PAHs, CO denotes ro- vibrational cool- 
ing by CO molecules, PE denotes photo-electric heating by PAHs, CR de- 
notes cosmic-ray heating, and PD denotes photo-dissociative heating of H2 
molecules. 

u cnv = 20 km s~\ comparable to the case where Hi medium 
is swept by a galactic spiral shock or by collisions of Hi shells 
in super-bubbles. We use periodic boundary condition at the 
y- and z-boundary planes. The initial setting is also periodic 
in the ^-direction except for v x given above. In effect, we im- 
pose continuous flows of the initial two-phase medium at the 
x -boundary planes at x = — 10 pc and 10 pc. Thus, there is 
a continuous mass input into the numerical domain from the 
x -boundaries. In this paper, we omit the dynamical effects of 
self-gravity, but we estimate it using a posteriori calculations 
in §3.2 and §3.5. The neglect of self-gravity for the duration of 
this simulation (t — 10.0 Myr) is justified, because the grav- 
itational energy does not become the dominant component of 
the energy budget (see §3.2) and most of the dense clumps 
in the molecular cloud stay in a gravitationally unbound state 
(see §3.5). 

Since the converging flow velocity is a few times larger 
than the sound speed in the WNM (c s w ~ 10 km s _1 ) and 
much larger than that in the CNM (~ 1 km s^ 1 ), two shock 
waves that bound the shocked slab are formed. In order to 
consider the effects of ultra-violet (UV) photon shielding and 
the resulting formation of molecules and the transition of 
the cooling/heating processes, we take into account the de- 
pendence of the local column densities (TVj) on the chemi- 
cal reaction coefficients and cooling/heating rates (see Tables 
1 and 2). In this paper, we assume that the sources of the 
UV field, which dissociates the molecules and causes photo- 
electric heating by PAHs, are located far from the region of 
molecular cloud formation. Because of the periodic condi- 
tions for ylz boundaries, the global morphology of the gener- 
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FIG. 2. — Volume rendering map of the density resulting from 8 Myr inte- 
gration (several cooling times for the initial thermally unstable gas). Regions 
in blue are Hi clouds with n > 10 cm -3 and T ~ 100 K that are formed 
as a consequence of the thermal instability. Magnetic field lines are shown 
as black lines. A two-dimensional slice is shown as the top surface in which 
red represents n < 3 cm -3 , yellow and green represent 3 cm~ 3 < n < 10 
cm— 3, and blue represents n > 10 cm -3 . The PDF of the number density 
is plotted in the bottom panel. 

ating molecular cloud should be regarded as sheet-like struc- 
ture that lies approximately in the y-z plane. Thus, the UV 
photons propagating along the y-z plane are well shielded in 
contrast to those propagating along the x direction. Hence, 
we use a two ray approximation for UV field that irradiates 
the numerical domain from the x = ±10 pc planes toward the 
=Fx directions with strength Go = 0.5 in both rays. Note that 
the shielding factor of the UV flux is larger for an isotropic 
background UV radiation field, since the shielding increases 
with increased column density in the pathway of UV pho- 
tons if the angle between a UV ray and the x-axis becomes 
larger. This means that the formation of molecules is eas- 
ier for an isotropic background UV radiation field with the 
same mean intensity. Thus, the two-ray approximation con- 
sidered in this analysis gives the upper bound for the for- 
mation timescale of molecules. We also calculate the ap- 
proximate escape probabilities for the emission of Cn/Oi fine 
structure lines and CO ro-vibrational lines using column den- 
sities along ±x directions, because the medium is opaque in 
the y- and z-directions. The local chemical reaction coeffi- 
cients (cooling/heating rates) are modified from the optically 
thin case by multiplying them by the shielding factors (escape 
probabilities) that are calculated by averaging the two shield- 
ing factors for the right and left rays (two escape probabilities 
for the left and right directions). Their specific forms, which 
were originally formulated in a slab geometry, are found in 
the references in Tables 1 and 2. 

In Figure [T|a), we plot the thermal and chemical equi- 
librium state under the assumed visual extinction of Av=l 
(iV tot = N K +2N Ka +N 1 ,+N He +N BB+ = 1.9xl0 21 cm- 2 ) 
as a dashed line. In panel (b) of Figure [T] we also plot the 
cooling and heating rates of the thermal and chemical equi- 
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librium states for the Av=0 (optically thin) and Av=l cases 
as solid and dashed lines, respectively. These panels shows 
that our numerical code can adequately treat media of warm 
atomic gas of T ~ 10 4 K to cold molecular gas of T ~ 10 K. 

3. RESULTS 

3.1. Formation of a Molecular Cloud 

The converging Hi flows induce two accretion shock waves 
that bound a highly inhomogeneous cloud. Owing to contin- 
uous mass accretion, the mean column density of the cloud 
along the x-axis increases linearly with time and a molecu- 
lar cloud is formed. In the top panel of Figure [3] we show 
the three-dimensional density map of the numerical domain 
at t = 10.0 Myr after the onset of the converging flows. Re- 
gions in yellow and green represent 3 cm~ 3 < n < 30 cm -3 
that correspond to the preshock Hi clouds and shocked WNM, 
regions in blue represent 10 2 cm~ 3 < n < 10 3 cm~ 3 that cor- 
respond to the shocked Hi clouds and molecular clumps, and 
regions in magenta represent dense clumps of n > 10 3 cm -3 
that are formed, e.g., by collisions of clumps. The PDFs of 
the number density (solid) and the initial HI medium (dotted) 
are plotted in the bottom panel. The PDF at t = 10 Myr is cal- 
culated only for gas in the region bound by the two accretion 
shocks. We show cross-sectional maps of the number density 
(n = Y^, n i) an d temperature in the z = and x = planes 
in Figure |4j where arrows represent the projected velocity on 
each plane. 

As shown by Koyama & Inutsuka (2002), even if the 
preshock medium is composed of WNM alone, nonlinear evo- 
lution of the thermal instability leads to the formation of tur- 
bulent cloud composed of small-scale cloudlets. In the present 
case, because the accreting two-phase medium is much more 
inhomogeneous than the single-phase WNM, the generated 
structure of the cloud is even more turbulent and inhomoge- 
neous. This is because the interaction of the accretion shock 
waves with upstream large density inhomogeneities triggers 
the Richtmyer-Meshkov instability (see, e.g., Nishihara et al. 
2010 for review) that induce rippling in the shock front. As a 
result, turbulent vorticity is created behind the curved shock 
as a consequence of Crocco's theorem (Kida & Orszag 1990; 
Kevlahan & Pudritz 2009; see also Inoue et al. 2009, 2010, 
2012 for recent simulations of the shock-cloud interaction). 
In the panel (b) of Figure [4] it is evident that the two accre- 
tion shock fronts have became curved due to the Richtmyer- 
Meshkov instability. Note that the growth timescale of the 
Richtmyer-Meshkov instability is given by the shock crossing 
time of the accreting CNM ~ 1 pc/20 km s _1 ~ 0.05 Myr, so 
it is influential from the very early stage of cloud formation. 

In the top panel of Figure[5] we plot the evolution of the av- 
erage column densities along the x-axis for H and He atoms 
(Nn = (J[n p + n H + 2n H2 + n Hc + rt He +]dx) ytZ ; solid), 
for H 2 molecules (N^ 2 = (J nn 2 dx) yiZ ; dashed), and for 
the CO-fraction-weighted density (Nqo = (J fco [ n p + 
Uh + 2jih 2 + «Ho + n He+]d x )y,z'-, dotted), where fco = 
nco/( n c+ + n c + n-co)- In the bottom panel, we plot the 
evolution of total mass (M to t = J P dV; solid), H 2 mass 
(Mh 2 = J 2 m p hh 2 dV; dashed), and CO-fraction-weighted 
total mass (Mco = J fcopdV; dotted) in the numerical 
domain. Because self-shielding of UV is effective for H 2 , 
whereas UV shielding by dust is necessary for CO formation, 
the H 2 mass and column density is always larger than the CO- 
fraction-weighted mass and column density. At t > 5 Myr, 
the total mass and the molecular mass measured by CO be- 
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FIG. 3. — Volume rendering map of the density at t = 10.0 Myr after 
the onset of converging flows (top). Regions in yellow and green represent 3 
cm~ 3 < n < 30 cm -3 that correspond to preshock Hi clouds and shocked 
WNM, regions in blue represent 10 2 cm -3 < n < 10 3 cm -3 that corre- 
spond to shocked Hi clouds and molecular clumps, and regions in magenta 
represent dense clumps of n > 10 3 cm -3 that are formed, e.g., by collisions 
of clumps. Magnetic field lines are shown as black lines. The PDF of the 
number density (solid) and that of the initial Hi medium (dotted) are plotted 
in the bottom panel. The PDF at t = 10 Myr is calculated only for gas in the 
region bounded by two accretion shocks. 




FIG. 4. — Panel (a): cross section map of number density in the z = 
plane. Panel (b): temperature map in the 2 = plane. Panel (c): number 
density map in the x = plane. Panel (d): temperature map in the x = 
plane. All the maps are snapshots at t = 10.0 Myr. Arrows represent the 
projected velocity on each plane. 

comes of the same order of magnitude. Thus, in the present 
case, the timescale at which most of C nuclei in the cloud 
become molecule is 5-10 Myr, which is consistent with the 
estimated timescale of molecular cloud formation associated 
with super-shells (Dawson et al. 2011). 

It has been pointed out that, in a cloud formed under the 
influence of the thermal instability, the cloud is not com- 
posed only of cold gas, but also the fragmented clumps of 
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FIG. 5. — Top: Average column densities along rc-axis for H and He 
atoms (JV H = (/[n p + n H + 2n Ha + riHe + 'n- ae +]dx) y y, solid), forH2 
molecules (7Vh 2 = (/ nn 2 dx) y , z ; dashed), and for CO-fraction-weighted 
density (N co = (/ fco [n p +n H +2 "H 2 + n Hc+n Hc +]dx) yiZ ; dotted), 
where fco = n CO I (p-c-v + + n Co)- Bottom: Evolution of total mass 
(Aftot = / pdV; solid), H2 mass (Mjj 2 = f 2in p nn 2 dV; dashed), and 
CO-fraction-weighted total mass (Aifco = / J CO P dV; dotted). 

clouds embedded in warm diffuse gas (T > 1000 K and 
n < 10 cm~ 3 ; e.g., Koyama & Inutsuka 2002, Hennebelle et 
al. 2008). Note that the diffuse warm gas is the thermally un- 
stable shocked WNM that eventually evolve into cold clumps 
due to the cooling (Banerjee et al. 2009). In previous studies 
of molecular cloud formation by converging WNM flows, the 
simulations are performed assuming optically thin radiative 
cooling and heating. Thus, the warm gas was irradiated by 
an overestimated UV field that led to excessive photoelectric 
heating and thus lengthened the cooling time and the lifetime 
of shocked warm gas. In the present case, however, even after 
the formation of molecules, i.e., even in the medium where the 
UV field is shielded, the shocked warm gas between the cold 
clumps survives as seen in Figure |4] The volume fraction of 
diffuse gas with T > 10 3 K ((n) T>10 3 K = 8.44 cirr 3 ) in the 
region bounded by the two accretion shock waves at t = 10.0 
Myr is 51.5%! The reason for the survival of this diffuse gas 
is explained as follows: The diffuse gas (more precisely the 
shocked WNM) is continuously supplied by the two accretion 
shocks. Since the typical cooling time for the shocked dif- 
fuse gas is only a few Myr, it seems very hard for the diffuse 
gas to exist deep inside the molecular cloud. However, as we 
will show in more depth in §3.3, the shock bounded region 
is highly turbulent and the velocity dispersion for the diffuse 
component of T > 10 3 K is Av > 10 km s -1 . Thus, the 
turbulent mixing length is evaluated as i m j x ~ Av t coo \ > 10 
pc, which makes the supply of diffuse gas to deep inside the 
molecular cloud possible before the gas is cools down. Note 
that, as discussed in Hennebelle & Inutsuka (2007), heating 
due to ambipolar diffusion makes the net cooling time longer 
than in the present case and would further increase the lifetime 
of diffuse gas in realistic molecular clouds (see also, Inoue et 
al. 2007 and Vasquez-Semadeni et al. 2011 for the effects of 
ambipolar diffusion on cloud formation). 

3.2. Evolution of Global Quantities 
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FIG. 6. — Evolution of the kinetic (red), magnetic (blue), thermal (green), 
and gravitational (magenta) energies. Kinetic and gravitational energies 
calculated using the CO-fraction-weighted density field are also plotted as 
dashed red and dashed magenta lines, respectively. The energies are calcu- 
lated only for the gas in the region bounded by the two accretion shock waves, 
i.e., energies in the unshocked converging flows are excluded. 

In Figure [6] we plot the evolution of the kinetic, ther- 
mal, magnetic, and gravitational energies as solid fines in 
red, green, blue, and magenta, respectively. We also plot 
the kinetic and gravitational energies using the CO-fraction- 
weighted density field as dashed red and dashed magenta 
lines, respectively. These energies are calculated only for 
gas in the region bounded by the two accretion shock waves, 
i.e., the energies in the unshocked converging flows are ex- 
cluded. The positions of the two accretion shocks (ir s h[j/, z}) 
are identified by searching for the maximum and minimum 
a; -coordinate values where the thermal pressure is larger than 
Pth = 10 4 fee er g cm~ 3 , since the converging Hi medium is 
almost isobaric with (p) = 5 x 10 3 fcs erg cm~ 3 . The gravi- 
tational energy (J p {r — r c } • V-0 dV) is calculated using the 
numerical solution of the Poisson equation (Aip = 4ir G p, 
and Ait G p fco f° r CO-fraction-weighted one) by assuming 
vacuum state outside the boundaries in x direction and pe- 
riodic boundary conditions for the y- and z-directions. The 
gravitational energies calculated for the raw density field and 
for the CO-fraction-weighted field become comparable (same 
order) at t ~ 7 Myr, which is consistent with the result shown 
in Figure [5] that the total mass and the CO-fraction-weighted 
mass become the same order at t ~ 5-10 Myr. The epochs 
at which the gravitational energy and turbulent energy mea- 
sured by the CO-fraction-weighted data become comparable 
are approximately the same (t ~ 7 Myr). This indicates that 
the scenario of molecular cloud formation examined in this 
paper is consistent with the fact that most observed molecular 
clouds are in an approximately virial equilibrium state. 

As discussed in the previous section, the resultant molec- 
ular cloud is composed of dense cold clumps embedded in 
diffuse warm gas. Figure [6] shows that the total kinetic en- 
ergy in the shock bounded region (solid red) is an order of 
magnitude larger than that contained in the molecular clumps 
(dotted red). This is because the velocity dispersion of the dif- 
fuse gas component is an order of magnitude larger than that 
of the molecular gas, while the density of the molecular gas 
is only an order of magnitude larger than the atomic gas and 
the volume filling fractions of the atomic and molecular com- 
ponents are comparable. For instance, in the shock bounded 
region at t = 10.0 Myr, the average density and velocity dis- 
persion in gas with CO fraction fco is larger than 0.1 are 
(n) / co >o.i = 383 cm~ 3 and Aw^ co>01 = 2.86 km s -1 , and 
in atomic gas with fco < 0.1 they are («)/ co <o.i = 35.0 
cm~ 3 and Au/ co> o.i = 12.2 km s _1 . This suggests that the 
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FIG. 7. — Evolution of the density-weighted (red), CO-density-weighted 
(green), and volume-averaged (blue) velocity dispersions in the shock 
bounded region. Solid lines represent the x component velocity dispersion 
(Av x = {{v x — (vx)} 2 ) 1 / 2 ), dashed lines represent that of the y and z 
component average ({Av 2 + Av 2 } 1 / 2 /2). A typical sound speed in molec- 
ular gas, c s ~ 0.3 (T/20K) 1 / 2 km s _1 , is also plotted as a dashed line, 
where the temperature T = 20 K is the average CO-density-weighted tem- 
perature. 

true turbulent kinetic energy in the cloud as a whole can be 
much larger than the kinetic energy of turbulence estimated 
using line -widths of molecular emissions. 

3.3. Turbulence 

In Figure [7] we plot the evolution of the density -weighted, 
CO-density-weighted, and volume averaged velocity disper- 
sions in the shock bounded region as red, blue, and green 
lines, respectively. Solid lines represent the m-component 
velocity dispersion (Av x = {(v x — (vx)} 2 ) 1 ^ 2 ), dashed 
lines represent those of the y- and z-component average 
({Av 2 +Av 2 } 1 ^ 2 /2). Atypical sound speed in molecular gas, 

c s ~ 0.3 (T/20 K) 1 / 2 km s"\ is also plotted as a dashed line, 
where the temperature T = 20 K is the CO-density-weighted 
average temperature at t = 10 Myr. In previous studies, it 
was shown that turbulence driven by the thermal instability 
(and the nonlinear thin shell instability) can induce "super- 
sonic turbulence" as a superposition of translative motions of 
clumps in the diffuse warm gas (Koyama & Inutsuka 2002), 
but the cold clumps have subsonic internal velocity disper- 
sions (Koyama & Inutsuka 2006, Heitsch et al. 2006). How- 
ever, in the present simulation, as we will show in §3.5.1, the 
velocity dispersion can be supersonic in cold clumps depend- 
ing on their size (see Figure |9l. This qualitative difference 
can be attributed to the additional generation of turbulence 
triggered by the interaction between the accreting CNM and 
the accretion shock, i.e., the Richtmyer-Meshkov instability, 
which was absent in previous studies of cloud formation. 

Figure [7] clearly shows that the turbulence in the molecu- 
lar cloud is anisotropic. This is quite reasonable because the 
orientation of the converging flows in the simulation is fixed 
that leads to the selective input of the source of the turbulence 
(see also Vazquez-Semadeni et al. 2007). It is widely known 
that most observed molecular clouds are filamentary, indicat- 
ing that the building blocks of molecular clouds are also ac- 
cumulated anisotropically. Thus, we can reasonably expect 
that turbulence in molecular clouds is anisotropic as well. Re- 
cently, Hansen et al. (201 1) reported that the dissipation rate 
of anisotropic turbulence can be smaller than that of isotropic 
turbulence, because the timescale of the cascade is determined 
by the crossing timescale in the direction of smaller velocity 
dispersion. As we have discussed in the previous sections, the 
cold clumps are embedded in diffuse warm gas with typical 



sound speed c s > 5 km s _1 , because the volume averaged 
temperature of the shock bounded region is (T) = 2280 K 
at t = 10.0 Myr. In isothermal gas, which is assumed in 
the conventional picture of molecular clouds, any supersonic 
motion of a cold clump creates a shock wave that deceler- 
ates the clump motion. However, since the sound speed in 
the warm gas is larger than the velocity of the cold clumps 
(~ velocity dispersion of the turbulence), the molecular cloud 
formed in this simulation is a system with less shock waves 
than the conventional isothermal model in which the dissi- 
pation rate of turbulence is diminished. Thus, owing to the 
effects of anisotropy and the diffuse warm gas, we can expect 
supersonic turbulence to be sustained longer, even after the 
accretion ceases. 

3.4. Magnetic Field 

In the present simulation, we set up accretion flows 
along the mean magnetic field, because a molecular 
cloud can only be formed when the raw materials are 
accumulated along the mean magnetic field (Inoue & 
Inutsuka 2009). This leads to a constant mean magnetic 
field strength (|(-B)| = 5.0 /iG) throughout the evolu- 
tion of the molecular cloud. Thus, the mass-to-flux ratio 
(fi = M/$ = 2 (p)v cnv t/\ (B) |) becomes larger than the 
critical value (|t cr = 0.13 G -1 / 2 ; Mouschovias & Spitzer 
1976) well before the formation of molecular cloud at t cr ~ 



2Myr((n)/5< 



(t/cnv^Okms- 1 )- 1 |/5 /xG). 



Note that the mass-to-flux ratio of the initial two-phase 
medium in the domain is [i//j, ci = 0.27. Even though 
the mean field is constant, magnetic field lines are locally 
stretched and compressed due to turbulence which results in 
a local amplification of the magnetic field. In panel (a) of 
Figure [8] we plot the average magnetic field strength ((\B\)) 
as a function of the local density at t = 5.0 Myr (red) and 
t = 10.0 Myr (blue), where bars indicate the dispersion of the 
magnetic field strength at that density. Clearly, the average 
magnetic field strength (\B\) in the dense gas is much larger 
than the mean magnetic field strength \(B)\ = 5.0 fiG, 
indicating that the magnetic energy of the molecular cloud 
shown in Figure [6] is dominated by the energy contained in 
magnetosonic waves rather than the energy of the global 
mean field. 

Because the turbulence in the formed molecular cloud is not 
only supersonic but also super-Alfvenic, it can easily bend 
and stretch the magnetic field lines. If the turbulence was 
isotropic, the passive nature of the magnetic field would re- 
sult in an isotropic distribution of local magnetic field orien- 
tations. However, as shown in the previous section, the tur- 
bulence is anisotropic and enhanced in the direction of the 
accretion flows. This leads to an anisotropic distribution of 
local magnetic field orientations. In panel (b) of Figure [8] 
we show the probability distribution function of the angle 
between the local direction of the magnetic field and the ini- 
tial direction (+x direction) in the shock bounded region at 
t = 10 Myr. The red line represents the density-weighted dis- 
tribution function, and lines in green, blue, magenta, and light 
blue represent the distribution functions in gas with 1 cm~ 3 < 

-3 < 



< n 

,-3 



< 10 2 



10 2 



10 4 cm 3 , respec- 



n < 10 cm -3 , 10 cm - ' 
n < 10 3 cm -3 , and 10 J cm~ J < n < 
tively. In panel (b), the plotted functions are normalized by 
dividing the calculated distribution by the isotropic distribu- 
tion Piso(9) oc sin(#), i.e., the flat distribution corresponds to 
the isotropic distribution. 
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FIG. 8. — Panel (a): Average strength of magnetic field as a func- 

tion of density at times t = 5.0 Myr (red) and t = fO.O Myr (blue). Bars 
indicate the dispersion of the magnetic field strength at a given density. Panel 
(b): Probability distribution function of the angle 9 between the local direc- 
tion of the magnetic field and the initial direction (+x direction) in the shock 
bounded region at t = 10 Myr. In panel (b), the plotted functions are nor- 
malized by dividing the calculated distribution by the isotropic distribution 
Piso(S) oc sin(f3), i.e., the flat distribution corresponds to an isotropic dis- 
tribution. 

The direction of the local magnetic field is biased towards 
the direction of the converging flows, except in high-density 
regions with n > 10 3 cm~ 3 . The degree of anisotropy de- 
creases with increasing density. This is because as the gas 
evolves towards increasing densities due to the condensation 
of thermally unstable gas and/or the collision of clumps, the 
local flow direction, which affects the orientation of the lo- 
cal magnetic field, is deflected from the initial direction of 
the converging flows. Note that the nearly isotropic distribu- 
tion in dense gas with n > 10 3 cm~ 3 does not indicate the 
random orientation of magnetic field in each dense clumps. 
It is noteworthy that, even in the low density regions with 
1 cm~ 3 < n < 10 cm~ 3 , the orientation of the magnetic field 
is flipped from the original direction (9 > 90 deg.) with cer- 
tain probability. The flipped orientation of the magnetic field 
even in the low-density gas reflects the fact that the postshock 
diffuse gas (shocked WNM) is not in a state of "laminar post- 
shock flow", but in a state of "anisotropic turbulence". 

3.5. Clump Statistics and Evolution 

In the formed cloud, molecular clumps evolve towards 
larger masses due to condensational accretion of the thermally 
unstable diffuse gas and collisional coalescence of clumps. In 
this section, we present the statistics and evolution of physical 
quantities of the clumps. In the following, we define a clump 
as a connected region with density is larger than a threshold 
value nth- In order to minimize the effects of numerical noise, 
we use the data for clumps that are resolved by more than 100 
numerical cells. 

3.5.1. Plasma f3, l-Sv relation, and mass-to-flux ratio 
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FIG. 9. — Scatter plots of the plasma f3 ((8 np/B 2 )) as attraction of clump 
mass and the velocity dispersion as a function of clump scale with threshold 
densities of n th = 100 cm -3 , 500 cm -3 , 2,000 cm -3 , and 10,000 cm -3 . 
Panels (a) and (b) show the plasma at t = 5.0 Myr and t = 10 Myr, 
respectively, and panels (c) and (d) show the size-velocity dispersion relation 
at t = 5.0 Myr and t = 10 Myr, respectively. The size of the clumps is 
evaluated as I = \J I\ /M, where I± is the largest eigenvalue of the inertial 
matrix and M is the mass of the clump. 

Figure|9]shows scatter plots of the plasma /3 (= (8ttp/B 2 )) 
as a function of clump mass, and plots of the velocity disper- 
sion as a function of clump scale. Points in red, green, blue, 
and magenta represent the clumps identified using thresh- 
old densities of n til = 100 cm~ 3 , 500 cirr 3 , 2,000 cm~ 3 , 
and 10,000 cm -3 , respectively. Panels (a) and (b) show the 
plasma f3 at t = 5.0 Myr and t = 10 Myr and panels (c) 
and (d) are the size- velocity dispersion relation at t = 5.0 
Myr and t = 10 Myr, respectively. The size of the clumps is 
evaluated as I = *J I\ JM, where Ii is the largest eigenvalue 
of the inertial matrix and M is the mass of the clump. We 
also used the geometric mean of the three eigenvalues and the 
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cubic root of the clump volume, but they did not change the 
results significantly. In panels (c) and (d), the typical sound 
speed in molecular gas, c s ~ 0.3 (T/20K) 1 / 2 km s _1 , is plot- 
ted as a thin dotted line, where the average temperatures in 
the clumps with n > 10 2 , 10 3 , and 10 4 cm -3 are 42, 19, 
and 13 K, respectively. The panels (a) and (b) indicate that 
the clumps evolve with time towards larger mass and that the 
number of high-density clumps increases with time, while the 
typical plasma f3 stays at (j3) ~ 0.4 almost independent of 
clump density and time. The size-velocity dispersion relation 
(panels [e] and [f]) also shows the universal law of 6v ~ 1.5 
km s _1 (7/1 pc) 5 , irrespective of both the mean density and 
time. This relation is compatible with the well known size- 
line width relation of molecular clouds (Larson 1981, Heyer 
& Brunt 2004), and consistent with the density independent 
relation observed in molecular cloud L1551 by Yoshida et al. 
(2010). 

3.5.2. Virial Parameters 

As clumps evolve, their mass, density, and size increase 
towards becoming a gravitationally unstable object. Panels 
(a) and (b) in Figure [10] show scatter plots of the clumps 
in a simple virial diagram in which the horizontal and ver- 
tical axes represent X = (U + K)/W and Y = M/W, 
respectively, where U = J 3pdV = 3 (7 — 1) E^, K = 
J P Sv 2 dV = 2E kia , M = J B 2 /8ndV = E mag , and 

W = J pr- Wip dV = E gvv , where the post-process gravita- 
tional potential ijj is used in the evaluation of the gravitational 
energy term W (see §3.2). Panels (a) and (b) show the results 
at t = 5.0 Myr and t = 10.0 Myr, respectively. A clump that 
enters the region of X + Y < 1 can evolve into star(s) due to 
gravitational contraction. Because clumps grow for f3 ~ 0.4 
or M ~ U irrespective of their density, the clumps lie approx- 
imately along the line X = Y for X > 1, For X ~ 1, the 
kinetic energy term K due to turbulence becomes comparable 
to or larger than the thermal energy term U. Thus, the clumps 
become gravitationally unstable cores with X > Y, indicat- 
ing that magnetically super-critical cores (whose thermal and 
turbulent kinetic energies are larger than the magnetic energy) 
are formed as a consequence of the clump evolution. 

Scatter plots of the mass-to-flux ratio ({E grv / E^g} 1 ^ 2 ~ 
ju//x cr ) as a function of clump mass at t = 5.0 Myr and 
t = 10.0 Myr are also shown in panels (c) and (d) of Fig- 
ure [TO] respectively. Clumps identified with the same n tn fol- 
low fi/fici oc M 0A , which is consistent with the result of 
Banerjee et al. (2009), while a higher nth leads to a larger 
mass-to-flux ratio for the same mass. When collisional co- 
alescence of clumps and accretion of thermally unstable gas 
onto the clumps occur along the magnetic field, the mass-to- 
flux ratio of clumps is enhanced. Typical high density clumps 
with n t h = 10 4 cm -3 can be super-critical when they grow 
to M > 1 M so i ar . It is widely known that the effects of 
ambipolar diffusion can enhance the mass-to-flux ratio even 
by condensational accretion of thermally unstable gas onto 
clumps perpendicular to the magnetic field (Inoue et al. 2007, 
Stone & Zweibel 2010) and when clumps are shocked by col- 
lisional coalescence (Li & Nakamura 2004, Kudoh & Basu 
2011). Thus, if we include the effects of ambipolar diffu- 
sion, the clumps can be super-critical even for M < 1 M so i ar . 
Note that clumps far from the unstable state X, Y ^> 1 (and 
M/Mcr "C 1) are small-scale, low-density, low -mass objects, 
and hence they should not be observed as dense magnetically 
subcritical cores with column density large enough to enable 
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FIG. 10. — Panels (a) and (b): Scatter plots of clumps in a simple virial 
diagram in which horizontal and vertical axes represent X = (U + K) /W 
and Y = M/W, respectively, where U = / 3pdV, K = / p6v 2 dV, 

M = J B 2 /8-KdV, and W = f pr- V0 dV. Panels (a) and (b) show the 
results at t = 5.0 Myr and t = 10 Myr, respectively. A clump that enters 
the region of X + Y < 1 (below the solid line) should evolve into star(s) 
by gravitational contraction. Panels (c) and (d): Scatter plots of the mass- 
to-flux ratio as a function of clump mass at t = 5.0 Myr and t = 10 Myr, 
respectively. 

Zeeman effect observations (e.g., Crutcher et al. 2009). 

3.5.3. Rotation of clumps 

Recent numerical simulations of the gravitational collapse 
of molecular cloud cores have revealed that cores can become 
fragmented after the firstcore is formed, if the initial condition 
of the core satisfies 



/p \l/2 > / (ffmag/ffth)^ 2 (\ / E mag / E th < 0.2) 

(iW-Vvj „ | Q 2 (^E mag /E th > 0.2) 

where E Iot , E mag , E t h, and E gl - V are the rotational, mag- 
netic, thermal, and gravitational energies of the clump 
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Panels (a) and (b) show the results at t = 5.0 Myr and t = 10.0 Myr, respec- 
tively. The dashed line represents the critical line of Machida et al. (2005) 
below which gravitational collapse leads to fragmentation into multiples. 
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solid line is the critical line for the above fragmentation con- 
dition. Panels (a) and (b) show the results at t = 5.0 Myr 
and t = 10.0 Myr respectively. Since even high-density 
threshold clumps with n t h = 10 4 cm~ 3 that are close to 
the state of gravitationally bound cores are located in the 
fragmentation region (above the dashed line), most clumps 
seem to evolve into cores that satisfy the fragmentation con- 
dition. The ratio of the rotational energy to the gravitational 
energy of observed molecular cloud cores is typically ~ 0.2 
or {E^/E^) 1 / 2 ~ 0.14 (e.g., Goodman et al. 1993) that 
suggests non-fragmentation for cores of -E mag /£ , th ~ 1. 
However, as shown in Dib et al. (2010), observations tend 
to underestimate the angular momentum by typically an or- 
der of magnitude which leads to comparable rotational en- 
ergies between the observations and our simulation. Note 
that, since neither the parameter E meig /Eth nor E rot /E gIV 
is a conserved variable, these parameters at the epoch of the 
beginning of the gravitational collapse are necessary to clar- 
ify whether the clumps truly evolve into cores that eventually 
fragment to form binaries, which requires a longer time sim- 
ulation including the effect of self-gravity. 

3.5.4. Mass spectrum 

Finally, we show the mass spectrum of the clumps. We plot 
the mass spectrum of the clumps at t = 5.0 Myr, and 10.0 
Myr in Figure 12 Panel (a) and (b) represent the spectra for 
n t h = 100 crrT^ and 1000 cm -3 , respectively. The mass 
spectrum of initial Hi clouds are plotted as dotted lines in 
both panels, where the threshold density is set to be nth = 20 
cm~ 3 (a few times smaller than the typical density of the ini- 
tial Hi clouds). The mass spectrum of the initial Hi clouds 
(dotted line) can be represented by a power law of the form 
dN/dM oc M' p with p ~ 1.8. This value of the spectral 
index p can be understood from the theoretical mass spec- 
trum of clumps formed by the thermal instability (Hennebelle 
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FIG. 12. — Mass spectrum of clumps at t = 5.0 Myr (solid), and 10 Myr 
(dashed). Panels (a) and (b) represent the spectra for n t h = 100 cm -3 and 
1,000 cm" 3 , respectively. The spectrum of initial Hi clouds is also plotted in 
both panels as a dotted line. 

& Audit 2007) in which the spectral index is calculated as 
p = 2 — (n — 3)/3 = 1.78, where n = 11/3 is the spectral 
index of the seed fluctuation grown by the thermal instabil- 
ity that is used in the initial condition of Hi cloud formation 
(see, §2.2). The spectra at later times show qualitatively dif- 
ferent aspects depending on the threshold density. The spec- 
trum with a lower threshold density ?i t h = 100 cm~ 3 (panel 
[a]) is similar to the spectrum of the initial Hi clouds, which 
is consistent with observations of molecular clumps (Kramer 
et al. 1996, Schneider 2002). At this lower threshold density, 
the spectra represent clumps generated by the compression of 
raw Hi clouds by the two accretion shocks. Thus, it is reason- 
able for the spectra in panel (a) to show similar shape to the 
initial Hi cloud spectrum. On the other hand, the higher den- 
sity threshold clumps in panel (b) show very different spectra 
from the initial one at both t = 5.0 and 10.0 Myr. Because 
the clumps identified using the high-density threshold would 
be more evolved, the deviation from the initial spectrum is not 
surprising. The shallower slope of the spectra for M < 1 M so i 
seems to suggest the growth of mass due to the coalescence 
of clumps by clump-clump collisions. It is noteworthy that, 
for M > 1 M so i, the spectra can be fitted by a power law with 
p = 2.3, which resembles the core mass function in molecular 
clouds (e.g., Ikeda et al. 2007, Ikeda & Kitamura 2009) and 
the high-mass tail of the initial mass function (Salpeter 1955), 
although the number of sample clumps is limited. 

4. SUMMARY AND DISCUSSION 

Using three-dimensional magnetohydrodynamic simula- 
tions including the effects of radiative cooling/heating, chem- 
ical reactions, and thermal conductivity, we have investigated 
the formation of molecular clouds in the ISM. We considered 
the formation of a molecular cloud due to the accretion of 
Hi clouds formed by the thermal instability. Such a mode of 
molecular cloud formation is compatible with the results of 



FORMATION OF TURBULENT AND MAGNETIZED MOLECULAR CLOUDS 



11 



recent observations (Blitz et al. 2006, Fukui et al. 2009). 
Since the mean density of the initial multiphase Hi medium 
is an order of magnitude larger than the typical diffuse WNM 
density, the formation timescale is shorter than in the previous 
studies of molecular cloud formation that assumed the accu- 
mulation of only diffuse WNM. The resulting timescale of 
molecular cloud formation of ~ 10 Myr is consistent with the 
evolutionary timescale of the molecular clouds in the LMC 
(Kawamura et al. 2009). 
The results of our simulation can be summarized as follows: 

• The initial multiphase Hi medium is compressed and 
piled up behind shock waves induced by accretion 
flows. Because the initial medium is highly inhomoge- 
neous, the accretion shocks are highly deformed. The 
deformed shock waves generate vortical motion behind 
the shock waves as a consequence of Crocco's theorem 
that results in the formation of a very turbulent molec- 
ular cloud. 

• The structure of the post shock region is roughly com- 
posed of dense cold gas (T < 100 K) and diffuse warm 
gas (T > 1, 000 K), and these two components are spa- 
tially well mixed by the turbulence (§3.1). Because the 
source of the turbulence is the accretion flows, the tur- 
bulence is highly anisotropic whose strength is biased 
toward the orientations of the accretion flows (§3.3). 

• The kinetic energy of the turbulence dominates the ther- 
mal, magnetic, and gravitational energies in the total 10 
Myr evolution. However, the kinetic energy measured 
by using the CO-fraction-weighted density is compa- 
rable to the other energies once the CO molecules are 
sufficiently formed owing to the UV shielding at 5-10 
Myr after the passage of shock wave, consistent with 
the fact that most of the molecular clouds are roughly 
virialized. This suggests that the true kinetic energy of 
the turbulence in the cloud as a hole can be much larger 
than the kinetic energy of the turbulence estimated by 
using line-width of molecular emissions (§3.2). 

• Since the turbulence is highly super-Alfvenic, local 
magnetic field strength in dense gas is amplified due 
to the so-called turbulent dynamo in addition to the 
compression, while the mean magnetic field strength 
in the cloud as a hole stays at the initial value. The 
local orientation of the magnetic field is also biased to- 
ward the direction of the accretion flows due to the ef- 



fect of anisotropic turbulence, except the dense gas with 
n > 10 3 cm" 3 (§3.4). 



• The clumps in the molecular cloud show statistically 
homogeneous evolution as follows: (i) The typical 
plasma (3 of the clumps is roughly constant ((3) ~ 0.4; 
(ii) The size-velocity dispersion relation show Aw ~ 
1.5 km s _1 (7/1 pc) 5 irrespective of density; (iii) the 
clumps evolve towards magnetically supercritical cores 
(see, Figure [T0| ); (iv) the clumps seem to evolve into 
gravitationally unstable cores that satisfy the fragmen- 
tation condition provided by Machida et al. (2005); (v) 
the mass function of the low density clumps can be fit- 
ted by a power law dN/dM oc M -18 that reflects the 
initial spectrum of the Hi clouds formed by the thermal 
instability, while that of the high density clumps sug- 
gests steeper slope of dN/dM oc M~ 2 3 for the high 
mass tail. 

In a subsequent paper, using synthetic observations of po- 
larized dust emission, we will show that the anisotropic orien- 
tation of magnetic field leads to an overestimation of the mean 
magnetic field strength when employing the Chandrasekhar- 
Fermi effect (Chandrasekhar & Fermi 1953). In this paper, 
we have omitted the dynamical effects of self-gravity. For the 
duration of our simulation {t = 10 Myr), the neglect of self- 
gravity can be justified, because the gravitational energy has 
not yet became the dominant component of the energy bud- 
get in the molecular cloud (see Figure [6]) and almost all dense 
clumps in the cloud are still in the gravitationally unbound 
state (see Figure 10 1. However, as clearly seen in Figure[6] the 
gravitational energy increases continuously, while other ener- 
gies plateau. Thus, in order to understand the later time evolu- 
tion, the inclusion of self-gravity is crucial (see Hennebelle et 
al. 2008, Heitsch et al. 2008, Banerjee et al. 2009, Vazquez- 
Semadeni et al. 201 1 for the cases of colliding WNM flows), 
which we will study in our next paper. 
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